#include "glm.hpp"
#include "ext.hpp"
#include<cstdio>
#include <iostream>
#include "PrintTem.h"
#include <gtx/string_cast.hpp>

//https://stackoverflow.com/questions/4124041/is-opengl-coordinate-system-left-handed-or-right-handed/12336360#12336360
//https://www.kevinbeason.com/smallpt/

using namespace glm;
#define Vec dvec3
struct Ray
{
    //远点 方向
    Vec o, d;
    Ray(Vec o_, Vec d_) : o(o_), d(d_) {}
};
enum Refl_t
{
    DIFF,
    SPEC,
    REFR
};
struct Sphere
{
    double rad;
    //位置 自发光 本色
    Vec p, e, c;
    Refl_t refl;
    Sphere(double rad_, Vec p_, Vec e_, Vec c_, Refl_t refl_) : rad(rad_), p(p_), e(e_), c(c_), refl(refl_) {}
    double intersect(const Ray &r) const
    {
        //ray origin --> sphere
        Vec op = p - r.o;
        double t, eps = 1e-4;
        double b = glm::dot(op, r.d);
        double det = b * b - glm::dot(op, op) + rad * rad;
        // print("intersect", glm::to_string(op), b, det, rad, glm::to_string(r.d));
        if (det < 0)
            return 0;
        else
            det = sqrt(det);
        return (t = b - det) > eps ? t : ((t = b + det) > eps ? t : 0);
    }
};
//radius position  emission color  material
// Sphere spheres[] = {
//     Sphere(1e5, Vec(1e5 + 1, 40.8, 81.6), Vec(), Vec(.75, .25, .25), DIFF),   //Left
//     Sphere(1e5, Vec(-1e5 + 99, 40.8, 81.6), Vec(), Vec(.25, .25, .75), DIFF), //Rght
//     Sphere(1e5, Vec(50, 40.8, 1e5), Vec(), Vec(.75, .75, .75), DIFF),         //Back
//     Sphere(1e5, Vec(50, 40.8, -1e5 + 170), Vec(), Vec(), DIFF),               //Frnt
//     Sphere(1e5, Vec(50, 1e5, 81.6), Vec(), Vec(.75, .75, .75), DIFF),         //Botm
//     Sphere(1e5, Vec(50, -1e5 + 81.6, 81.6), Vec(), Vec(.75, .75, .75), DIFF), //Top
//     Sphere(16.5, Vec(27, 16.5, 47), Vec(), Vec(1, 1, 1) * 0.999, SPEC),       //Mirr
//     Sphere(16.5, Vec(73, 16.5, 78), Vec(), Vec(1, 1, 1) * .999, REFR),        //Glas
//     Sphere(600, Vec(50, 681.6 - .27, 81.6), Vec(12, 12, 12), Vec(), DIFF)    //Light
// };

//HDR 亮度12的灯
Sphere spheres[] = {
    // Sphere(16.5, Vec(0, 0, 0), Vec(12, 12, 12), Vec(), DIFF),       //单Light

    //暗部更多颜色
    Sphere(26, Vec(0, 26+18, 0), Vec(12, 12, 12), Vec(), DIFF),       //单Light
    // // Sphere(5, Vec(20, 0, 0), Vec(0, 0, 0), Vec(0.75, 0.25, 0.25), DIFF),       //单Light

    Sphere(1e5, Vec(-1e5-20, 0, 0), Vec(0, 0, 0), Vec(0.75, 0.25, 0.25), DIFF),       
    Sphere(1e5, Vec(1e5+20, 0, 0), Vec(0, 0, 0), Vec(0.25, 0.75, 0.75), DIFF),       
    Sphere(1e5, Vec(0,-1e5-20, 0), Vec(0, 0, 0), Vec(0.25, 0.75, 0.25), DIFF),       
    Sphere(1e5, Vec(0,1e5+20, 0), Vec(0, 0, 0), Vec(0.75, 0.75, 0.75), DIFF),       
    Sphere(1e5, Vec(0, 0, -1e5-20), Vec(0, 0, 0), Vec(0.75, 0.75, 0.75), DIFF),       
    Sphere(5, Vec(-10, -15, 0), Vec(0, 0, 0), Vec(1.0, 1.0, 1.0)*0.999, SPEC),       
    Sphere(8, Vec(10, -15, 0), Vec(0, 0, 0), Vec(1.0, 1.0, 1.0)*0.999, REFR),       

    // Sphere(8, Vec(0, 0, 0), Vec(0, 0, 0), Vec(1.0, 1.0, 1.0)*0.999, REFR),       
};

double clamp(double x) { return x < 0 ? 0 : x > 1 ? 1 : x; }
//gamma 修正线性 到 gamma 空间存储颜色 向上取整
int toInt(double x) { return int(pow(clamp(x), 1 / 2.2) * 255 + .5); }
//Ray 相交 某个Sphere
bool intersect(const Ray &r, double &t, int &id)
{
    double n = sizeof(spheres) / sizeof(Sphere);
    double d; 
    auto inf = t = 1e20;
    for (int i = int(n); i--;)
    {
        if ((d = spheres[i].intersect(r)) && d < t)
        {
            t = d;
            id = i;
        }
    }
    return t < inf;
}

//相交某个Sphere 计算辐射度
/*
one ray everything 
one ray pixel
one shadow ray  for point light
one reflection ray  one refraction ray  per intersection
*/
Vec radiance(const Ray &r, int depth)
{
    double t;   // distance to intersection
    int id = 0; // id of intersected object
    // print("radiance", glm::to_string(r.o), glm::to_string(r.d), depth);
    if (!intersect(r, t, id))
        return Vec();                // if miss, return black
    const Sphere &obj = spheres[id]; // the hit object
    //相交点 法向
    Vec x = r.o + r.d * t; //交点
    //球体中心到交点方向 法向
    Vec n = normalize(x - obj.p); //交点处法向
    //修正法向和ray 朝向  法向量方向 和 r.d 锐角  不可能打出锐角 球体相交法向 总是和ray钝角 打不到背面 总是normal 反向
    //正确法方向 光线相交表面法线 < 0
    Vec nl = dot(n, r.d) < 0 ? n : n * -1.0;
    //表面颜色
    Vec f = obj.c;
    //    double p = f.x>f.y && f.x>f.z ? f.x : f.y>f.z ? f.y : f.z; // max refl
    //最大颜色 深度>5 随机Xi < p f = f*(1/p)  自发光 自己颜色最强的除法贡献
    //随机数
    //    if (++depth>5) if (erand48(Xi)<p) f=f*(1/p); else return obj.e; //R.R.
    //深度太深了 停止pathTrace
    //https://computergraphics.stackexchange.com/questions/2316/is-russian-roulette-really-the-answer
    if (++depth > 5)
    {
        auto rd = glm::linearRand(0.0, 1.0);
        double p = glm::max(f.x, f.y, f.z);
        if (rd < p)
        {
            f = f * (1 / p); //本原色标准化 后面使用 概率扩大的能量
        }
        else
        {
            return obj.e; //自发光颜色 丢失的能量
        }
    }

    //普通diffuse材质
    if (obj.refl == DIFF)
    { // Ideal DIFFUSE reflection
        //2pi *Xi 随机
        //随机散射方向 竟然是多次 uint8 拼起来的 uint64 位整数
        auto rd = glm::linearRand(0.0, 2 * M_PI);
        auto rd2 = glm::linearRand(0.0, 1.0);
        //0 ~ 2PI 随机
        double r1 = rd;
        //0~1 随机
        double r2 = rd2;
        //开根 0~1随机
        double r2s = sqrt(r2);
        //法向
        Vec w = nl;
        //up 方向 y轴 或者z轴  normal > 0.1x
        Vec u = fabs(w.x) > .1 ? Vec(0, 1, 0) : Vec(1, 0, 0);
        //cross Z 方向
        u = normalize(cross(u, w));
        //w u v  Cross Y方向
        Vec v = cross(w, u);
        //up 分量 + forward 分量 + normal 分量
        //表面坐标轴系统随机反射方向 w 相反半球面 随机一个反射光线   u, v, w 做个坐标轴 theta phi
        //(1-r2)^2 + r2s^2*cos(r1)^2 + r2s^2*sin(r1)^2 = 1
        // (u, v, w) 坐标轴  * u v 平面 2pi 角度 w 高度 0~1
        Vec d = normalize(u * cos(r1) * r2s + v * sin(r1) * r2s + w * sqrt(1 - r2));
        //  return obj.e + f.mult(radiance(Ray(x,d),depth,Xi));
        //交点位置 随机一个反射方向
        //自发光 + 随机吸收自身色 * 入射光线 各分量调整颜色 入射光线颜色 * 自身颜色
        return obj.e + f * radiance(Ray(x, d), depth);
    }
    else if (obj.refl == SPEC)
    { //理想镜面 完全反射          // Ideal SPECULAR reflection
        //  return obj.e + f.mult(radiance(Ray(x,r.d-n*2*n.dot(r.d)),depth,Xi));
        // return obj.e + dot(f, radiance(Ray(x, r.d-n*2*n)));
        //反射方向  cos normal direction  r+ 2* normal 投影
        auto relDir = r.d - 2.0 * dot(n, r.d) * n;
        auto relRay = Ray(x, relDir);
        //specular color 受 本原色影响
        return obj.e + f * radiance(relRay, depth);
    }
    //玻璃球
    auto reflDir = glm::reflect(r.d, n);
    Ray reflRay(x, reflDir);
    //    Ray reflRay(x, r.d-n*2*n.dot(r.d));     // Ideal dielectric REFRACTION
    //    bool into = n.dot(nl)>0;                // Ray from outside going in?
    auto into = glm::dot(n, nl) > 0; //入射光线
                                     //真空折射率
    double nc = 1;
    //玻璃折射率
    double nt = 1.5;
    //入射光线 折射率? 1/1.5 空气到玻璃 玻璃glass 折射率
    double nnt = into ? nc / nt : nt / nc;
    //光线法线投影 
    double ddn = glm::dot(r.d, nl); 
    //入射光线
    //1-折射率平方*(1-法线投影)
    //透射比 低于一定角度 完全反射 投射 反射比例
    double cos2t = 1 - nnt * nnt * (1 - ddn * ddn);
    //玻璃进入空气
    //nc/nt = sin(r)/sin(i)
    //sin(r) = sin(i) * nc/nt
    //sin(r) == 1 90度 无折射
    //1 - nnt * sin(i) < 0 无折射
    //sin(i)^2 = 1-cos(i)^2 = 1- dot(d, nl)^2
    //1 - sin(i)^2*nnt^2  全反射
    if (cos2t < 0)
    { // Total internal reflection
        //  return obj.e + f.mult(radiance(reflRay,depth,Xi));
        return obj.e + f * radiance(reflRay, depth);
    }
    //反射+透射
    //透射方向 normal + short
    //dot r.d, -nl
    // auto inLen = glm::dot(r.d, -nl);
    // auto chuiZhiDir = r.d - inLen * -nl;
    // auto shortDir = chuiZhiDir * nnt;
    // auto lenShort = glm::dot(shortDir, shortDir);
    // auto normalShort = (1 - lenShort);
    // auto refracDir = normalShort * -nl + shortDir;
    // //透射方向
    // auto tdir = glm::normalize(refracDir);
    auto tdir = glm::normalize(glm::refract(r.d, nl, nnt));
    //    Vec tdir = (r.d*nnt - n*((into?1:-1)*(ddn*nnt+sqrt(cos2t)))).norm();
    //透射比 R
    //https://en.wikipedia.org/wiki/Fresnel_equations#Power_(intensity)_reflection_and_transmission_coefficients
    //入射角 出射角 0度 cos theta = 1 近似反射能量
    auto a = nt - nc;
    auto b = nt + nc;
    //反射能量Rs
    auto R0 = a * a / (b * b);
    //入射 反向投影 出射 投射方向 在 法线上投影
    //cos theta = 入射角 cos theta
    //specular reflection coefficient R  反射能量 玻璃内能量反射
    auto c = 1 - (into ? -ddn : glm::dot(tdir, n));
    //    auto c = 1-(into?-ddn:tdir.dot(n));
    //https://en.wikipedia.org/wiki/Schlick%27s_approximation
    //R0 = (a/b)^2
    //R(theta) = R0 + (1-R0)(1-cos(theta))^5
    //反射能量Rp
    auto c5 = glm::pow(c, 5);
    auto Re = R0 + (1 - R0) * c5;
    //    auto Re=R0+(1-R0)*c*c*c*c*c;
    //透射能量 Tp
    auto Tr = 1 - Re;
    //过深
    if (depth > 2)
    {
        //概率 只考虑反射或者折射 概率 根据反射自身权重决定
        //反射强度概率
        auto P = .25 + .5 * Re;
        //权重/概率 修正权重
        auto RP = Re / P;
        //修正 透射权重
        auto TP = Tr / (1 - P);
        auto rdd = glm::linearRand(0.0, 1.0);
        //只考虑反射
        if (rdd < P)
        {
            return obj.e + f*radiance(reflRay, depth) * RP; //只考虑反射能量
        }
        else
        {
            return obj.e + f*radiance(Ray(x, tdir), depth) * TP; //只考虑透射
        }
    }
    else
    { //普通
        //反射+投射能量 比例
        return obj.e + f * (radiance(reflRay, depth) * Re + radiance(Ray(x, tdir), depth) * Tr);
    }
    //    return obj.e + f.mult(depth>2 ? (erand48(Xi)<P ?   // Russian roulette
    //      radiance(reflRay,depth,Xi)*RP:radiance(Ray(x,tdir),depth,Xi)*TP) :
    //      radiance(reflRay,depth,Xi)*Re+radiance(Ray(x,tdir),depth,Xi)*Tr);
}
void Test(){
    Vec camCen(50, 52, 295.6);
    Vec sphereCen(50, 16.5, 47);
    auto dir = sphereCen-camCen;
    auto d = glm::normalize(dir);
    // print(glm::to_string(camCen), glm::to_string(sphereCen), glm::to_string(d));
}
//测试样例 输入 不同的输出
//小规模 到大规模
/*
测试: 单个光源 分辨率逐渐增加
增加随机性路径
*/
int main()
{
    Test();
    //输出像素
    // auto w = 1024;
    // auto h = 768;
    // auto w = 10;
    // auto h = 10;
    // auto w = 1024;
    // auto h = 768;
    // //每个像素采样ray
    // auto samps = 2;
    // auto subPixel = 2;

    auto w = 1024;
    auto h = 768;
    // auto w = 3;
    // auto h = 3;
    //每个像素采样ray
    auto samps = 2;
    auto subPixel = 2;

    //镜头位置和方向
    // Ray cam(Vec(50, 52, 295.6), glm::normalize(Vec(0, -0.042612, -1)));
    Ray cam(Vec(0, 0, 200), glm::normalize(Vec(0, 0, -1)));
    //Near  200 0  100位置 -30 30 范围的镜头
    
    //宽度一半/h
    //摄像机空间-1 1 -1 1范围-->0 w范围
    // Vec cx(w * 0.5135 / h, 0, 0);
    //镜头前1m 距离的对应的世界高度
    auto fullScreenHeight = 0.3;// 1/200 = x/30 x = 0.15  2*x = h
    //整个屏幕的1距离的左右宽度和高度
    Vec cx(w*fullScreenHeight/h, 0, 0);

    //镜头方向 cx x方向  cy 屏幕Y 方向
    // Vec cy = glm::normalize(glm::cross(cx, cam.d)) * 0.5135;
    //整个屏幕的1距离高度
    Vec cy = glm::normalize(glm::cross(cx, cam.d)) * fullScreenHeight; //Y 方向每个像素对应空间距离
    print("cx cy", glm::to_string(cx), glm::to_string(cy));
    //放大距离
    auto screenPos = 100.0;
    auto subRate = 1.0/(subPixel*subPixel);

    Vec r;
    //输出像素
    Vec *c = new Vec[w * h];
    //逐行渲染
    //左下渲染 PPM 左上0

    #pragma omp parallel for schedule(dynamic, 1) private(r)
    for (int y = 0; y < h; y++)
    {
        printf("RenderLine %d %d %f\n", y, samps * subPixel*subPixel * y * w,  100.0f* y/(h));
        //每列
        for (int x = 0; x < w; x++)
        {
            //2*2
            //i 像素编号
            for (int sy = 0, i = (h - y - 1) * w + x; sy < subPixel; sy++)
            {
                //收集每个像素 2*2子像素和采样次数的结果
                for (int sx = 0; sx < subPixel; sx++, r = Vec())
                {
                    //每个子像素采样次数
                    for (int s = 0; s < samps; s++)
                    {
                        double r1 = glm::linearRand(0.0, 2.0);
                        auto r2 = glm::linearRand(0.0, 2.0);
                        //0.0~2.0 0.0 1.414 
                        //-1~0 0~1 但是分布式 sqrt 分布
                        //修改了随机分布模型 -1~1 sqrt 概率分布
                        auto dx = r1 < 1 ? sqrt(r1) - 1 : 1 - sqrt(2 - r1);
                        auto dy = r2 < 1 ? sqrt(r2) - 1 : 1 - sqrt(2 - r2);
                        //x基础像素位置 + sx 偏移 2*2子像素 +0.5 中心 + dx 随机采样 概率分布
                        //sx 0~1 dx -1~1 -1~2  w总宽度
                        //2*2 1*1 0.5 0.5 -0.25 + 0.75  0.25-0.5 
                        //范围-0.5 ~ 0.5 屏幕范围 对应的空间范围
                        // auto xmod = ((sx+0.5+dx)/2+x)/w-0.5;
                        // auto ymod = ((sy+0.5+dy)/2+y)/h-0.5;
                        // auto xmod = x*1.0/w -0.5 + 0.5/w; //子像素中心
                        // auto ymod = y*1.0/h - 0.5+0.5/h;
                        auto xmod = x*1.0/w -0.5 + (sx+0.5+dx/2.0)/(subPixel*w); //超采样2*2
                        auto ymod = y*1.0/h - 0.5+ (sy+0.5+dy/2.0)/(subPixel*h); 
                        //规范化空间 镜头前-1 位置
                        // xmod = x*1.0/w-0.5+0.25; //-0.25 0.25
                        // ymod = y*1.0/h-0.5+0.25; //-0.25 0.25
                        auto d = xmod * cx + ymod*cy + cam.d;
                        // print("direct:##", i, glm::to_string(d), xmod, ymod, glm::to_string(cx), glm::to_string(cy));
                        //采样次数平均
                        r = r+radiance(Ray(cam.o+d*screenPos, glm::normalize(d)), 0) *(1.0/samps);
                        //像素 对应的2*2 子像素 + 随机收集信息 
                        // Vec d = cx * (((sx + .5 + dx) / 2 + x) / w - .5) +
                        //         cy * (((sy + .5 + dy) / 2 + y) / h - .5) + cam.d;
                    }
                    // print("ReadColor:", glm::to_string(c[i]), glm::to_string(r));
                    c[i] = c[i] + Vec(clamp(r.x), clamp(r.y), clamp(r.z))*subRate;
                }
            }
        }
    }
    //输出PPM格式
    FILE *f = fopen("image.ppm", "w");
    fprintf(f, "P3\n%d %d\n%d\n", w, h, 255);
    auto cv = 0;
    for(auto i = 0; i < w*h; i++){
        fprintf(f, " %3d %3d %3d  ", toInt(c[i].x), toInt(c[i].y), toInt(c[i].z));
        cv++;
        if((cv % w) ==0){
            fprintf(f, "\n");
        }
    }
    fclose(f);
}